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Non-adiabatic effects play an important role in many chemical processes. In order to study the 
underlying non-adiabatic potential-energy surfaces (PESs), we present a locally-constrained density- 
functional theory approach, which enables us to confine electrons to sub-spaces of the Hilbert space, 
e.g. to selected atoms or groups of atoms. This allows to calculate non-adiabatic PESs for defined 
charge and spin states of the chosen subsystems. The capability of the method is demonstrated by 
calculating non-adiabatic PESs for the scattering of a sodium and a chlorine atom, for the interaction 
of a chlorine molecule with a small metal cluster, and for the dissociation of an oxygen molecule at 
the Al(lll) surface. 
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I. INTRODUCTION 

Assuming that electrons can react much faster to ex- 
ternal perturbations than nuclei, the Born-Oppenheimer 
approximation 1 (BOA) is a widely employed approach to 
separate electronic and nuclear motion in dynamic pro- 
cesses. The nuclei then appear static to the electrons, 
which in turn set up a potential governing the motion of 
the nuclei. Formally the approach starts with the general 
many-body Schrodinger equation 

HV = (T nuc + V nuc - nuc + JT e )¥ = EV , (1) 

where T nuc is the kinetic energy operator of the nuclei, 
ynuc-nuc ^ e interaction potential of the nuclei, and H e 
the electronic Hamiltonian containing the electron ki- 
netic energy, as well as the electron-electron and electron- 
nuclei interactions. In an adiabatic representational^, 
the general dependence of the full many-body wavefunc- 
tion ^({ridi}, {R/<7f }) on the position and spin coordi- 
nates of all electrons, r.; and o^, and on the position and 
spin coordinates of all nuclei, Rj and 07, is written as 

* = £ A^ dia ({R /( r J })^ ia ({r,c7 i },{R /< 7 / }) , (2) 

where the $^ dla are chosen to be the eigenfunctions of 
the electronic Hamiltonian at the actual position of the 
nuclei 

^e^adia = £e $ adia (3) 

Inserting Eq. ([2]) into Eq. ([I]) leads to a set of equations 
for the wavefunctions of the nuclei 

^^adia ^ynuc _|_ ynuc— nuc _j_ ^adia 

+non-adiabaticity terms , (4) 

in which the "non-adiabaticity terms" summarize matrix 
elements of the momentum and kinetic energy operators 



of the nuclear motion. The Born-Oppenheimer approxi- 
mation corresponds to setting these terms to zero, which 
implies that the electrons assume their electronic state in- 
stantaneously for any position of the nuclei, unaffected by 
the nuclear dynamics. The nuclei are then moving on the 
potential-energy surface (PES) V* dia = V nuc ~ nuc + Ef y7 
also called the Born-Oppenheimer surface of the elec- 
tronic state v. 

If, in the BOA, the system is initially in the electronic 
ground state, it will remain there irrespective of the dy- 
namics of the nuclei. However, in real life electrons may 
not be able to follow the motion of the nuclei instanta- 
neously, and, e.g. when selection rules apply, they may 
find themselves in an excited state. For example, chem- 
ical reactions forming singlet molecules from triplet and 
singlet reactants are forbidden by Wigner's spin selec- 
tion rule. And the triplet multiplicity is the actual rea- 
son, why most reactions of O2 with other molecules or 
substances, although being exothermic, do not proceed 
at room temperature; they are kinetically hindered. In 
other words, in a chemical reaction the spin of the reac- 
tants must be conserved or transferred to some other en- 
tity. And the transition from the O2 ground state, ( 3 ^~) 
to the first excited state ( 1 A S ), is strictly forbidden for 
the isolated molecule, as is the reverse (de-excitation) 
process, once the molecule has been excited to the X A 9 
state. Thus, when probabilities for transitions between 
different electronic states are low, e.g. due to selection 
rules, the assumption that the system will always remain 
in the electronic ground state may become incorrect. For 
the mentioned O2 example this implies that when an ex- 
ternal field shifts the 1 A S energy below the 3 S~ energy, 
the probability for an electronic transition toward the 
1 A g state will be low. Indeed, this is an important as- 
pect of the 02/Al(lll) interaction, one of our examples 
discussed below. 

For such, or alike, situations it is necessary to go be- 
yond the BOA and to consider the coupling between dif- 
ferent states££ For a full description the non-adiabaticity 
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terms must be calculated, which in practice means to 
evaluate the matrix elements for the derivative couplings 
of the nuclear and electronic motion. For this, the use 
of another set of electronic functions |$^ la > in the wave 
function expansion in Eq. ([2J) may be suitable. The idea 
behind such "diabatic" states is that in a dynamic (scat- 
tering) event, the electrons lag behind the nuclear motion 
and thus tend to conserve a given (initial) character of 
the electronic structure. Formally, the |$^ ia > are then 
constructed to minimize the nuclear derivative coupling 
terms, such that the matrix representation of the nuclear 
momentum operator becomes diagonal in the diabatic 
basis&££. The equivalent to Eq. ([4| in a diabatic repre- 
sentation is thus a matrix equation 

= E,A^ . (5) 
where the diagonal elements of the potential, V^ m — 

ynuc-nuc + jje^ &re ^ diabatic p ES g ; and the off _ 

diagonal elements, H^ v , describe the coupling between 
PES [a and PES v. For the above example of an O2 
molecule, a suitable diabatic state would e.g. be given 
by the 3 £~ state of the isolated molecule, and the corre- 
sponding diabatic PES can be used to describe the mo- 
tion of a molecule that tries to conserve this electronic 
configuration, even in regions where external fields (or 
interactions with other species) bring the 1 H g below the 
3 E~ energy. 

The aim of this work is to develop and apply a non- 
adiabatic approach, with a similar physical motivation 
as the diabaticity concept. We account for the tendency 
to maintain a given character of the electronic structure, 
e.g. due to selection rules, by focusing on the dynam- 
ics of a system with defined charge and spin in suitably 
chosen subsystems, typically the two reactants in a scat- 
tering event. Using again the example of an O2 molecule, 
this could e.g. be a state with fixed triplet spin on the 
molecule, even when it interacts with another reactant. 
Instead of using properties of the many-body wavefunc- 
tion as formal basis, we thus build our approach on a 
defined total spin (or total charge, or other well defined 
quantity) in the local Hilbert space of a chosen subsys- 
tem, which can be an atom or a group of atoms. The thus 
defined non-adiabatic PESs are computed with density- 
functional theory (DFT)^ 1 ^, and in particular using the 
idea of constrained DFT formulated by Dederichs et al. 
in 1984ii. The basic idea is to modify the energy func- 
tional employed in DFT by applying physically mean- 
ingful constraints. The constraint is enforced through 
an additional Lagrange multiplier, which on the level of 
the Kohn-Sham (KS) equations yields an additional con- 
straint term to the effective KS potential. In the present 
work we employ a local constraint via a projection tech- 
nique that enables us to freeze the charge and spin states 
of the chosen subsystems. 

Having stated this general philosophy, two points de- 
serve a critical discussion. First, we point out that in 



practice our non-adiabatic states may often be quite close 
to those of the diabaticity concept. Still, the formal def- 
inition is different, and in some cases like the example 
of the 02/Al(lll) system discussed below, there are no- 
table differences. In contrast to recent other works in 
this are a 12 ' 13 , we therefore prefer to refer to our states 
simply as "non-adiabatic states" to underscore this differ- 
ence to the established diabaticity concept. As a second 
point, we note that a solid, mathematical proof of the 
validity of constrained DFT does not exist. However, 
spin-density functional theory, the Slater- Janak transi- 
tion state concep t 11 ! 14 ' 15 , the fixed-spin moment (FSM) 
approach 1 » and some other examples present frequently 
employed, important, and successful applications of the 
method. Our concept represents a mild generalization 
of such applications, and our constraint is plausible and 
physically meaningful: The spins or charges of two in- 
teracting systems may be hindered to adjust or combine 
when there are selection rules. Thus, the appropriate 
total energy surfaces of the scattering event should (up 
to a certain distance) be that of conserved local spins or 
charges. 

In a previous publication we have already used the 
present approach to resolve a long-standing problem in 
surface science, namely the low sticking probability of 
oxygen molecules at the Al(lll) surfaced However, the 
method is much more general, as we will illustrate below 
by also applying it to two other, non-periodic model sys- 
tems, namely a scattering event of a sodium and a chlo- 
rine atom, and the interaction of a CI2 molecule with 
a magnesium cluster. In this respect, we also mention 
notable, independent work of Wu and Van Voorhi o 13 ' 18 , 
which is essentially equivalent to our approach^ in that 
also in their work an additional potential is introduced 
to constrain electron numbers in well-defined parts of a 
system. As in our approach this potential is determined 
in a self-consistent way using formally identical method- 
ology, and has been employed to study charge-transfer 
reactions. 



II. LOCALLY-CONSTRAINED 
DENSITY-FUNCTIONAL THEORY 

Our locally-constrained DFT (LC-DFT) approach 
starts by assigning electrons to defined sub-spaces of the 
total Hilbert space, e.g. to atoms or groups of atoms. 
This is done by employing a suitable projection scheme 
to distinguish the individual subsystems. In the follow- 
ing sections we will present this formalism for a sys- 
tem consisting of two subsystems called A and B, with 
straightforward generalization to more than two subsys- 
tems. Considering electrons and their spin, this leads 
to four electron numbers N^, 2Vj[, iVg and N^, which 
uniquely define the PES U™.^ ({ R '^}) for the 
non-adiabatic quantum state with fixed spins and charges 
of the subsystems. Expressing the constraints in terms 
of auxiliary potentials, the electronic structure problem 
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is solved self-consistently using DFT, i.e., the electronic 
structure is fully relaxed under the given constraints 
yielding Vfi^^. 

A. Definition of the Subsystems 

In order to assign electrons to the two subsystems, the 
Kohn-Sham single-particle wavefunctions are expanded 
into localized, atom-centered basis functions, e.g. Gaus- 
sians, Slater type orbitals, or numerical orbitals. The 
implementation is therefore particularly convenient in 
codes employing such basis sets, whereas in codes based 
on other basis sets, like plane waves, the implementa- 
tion requires intermediate projection steps onto localized 
functionsi&^^k The Kohn-Sham orbital <pf with index 
i and spin index a = f , { is thus written as a linear com- 
bination of all basis functions \j j or m calculations using 
periodic boundary conditions as a linear combination of 
k-dependent Bloch basis functions = e lkr Xj, 

n 

ft = £#tf ■ (6) 

3=1 

In the following our derivation will refer to the periodic 
case, while for finite systems the dependence on the k- 
point index would be simply dropped. 

In the case of atom-centered basis functions each basis 
function is uniquely assigned either to subsystem A or to 
subsystem B: All basis functions centered at atoms being 
part of subsystem A define the Hilbert space of subsystem 
A, and all basis functions centered at atoms being part 
of subsystem B define the Hilbert space of subsystem B. 
Every single-particle wavefunction can then be projected 
onto the two Hilbert sub-spaces, which is done separately 
for each k-point and each spin, and taking into account 
a non-orthogonality of the atomic orbitals by including 
the overlap matrix elements S^ k =< Xj \x\ >) 

n rn 

pZ = < ft\ft A > = EE^W ( 7a ) 

3=1 k=l 
n n 

Pti = < > = E E <%Sjk& • ( 7b ) 

j— 1 fe=m+l 

Here, p^ i is the projection onto subsystem A, p^f i the 
projection onto subsystem B, and n is the total num- 
ber of basis functions, of which the first k = 1 . . . m 
are the basis functions of subsystem A. <j^" A is defined 
as the A-component of <p\ a , i.e., all coefficients refer- 
ring to basis functions of subsystem B are set to zero. 
Correspondingly, all coefficients referring to basis func- 
tions of subsystem A are set to zero in the functions cf)^ 
with i = m + 1 . . . n, which define the complementary 
B-component of <j%" . We note that the normalization 
condition p^ i + p l ^ r i = 1 holds for each state i by con- 
struction, because each 4>f^ + 4>^% = <t>\ a is a normalized 



eigenstate. This relation can be used to reduce the com- 
putational effort in that just one of the two double sums 
in Eq. ((?]) needs to be calculated and the other is ob- 
tained from the normalization condition. We note that 
in this scheme like in all projection schemes the actual 
values of p Ai and p^ r i depend on the choice of the basis 
sets defining the Hilbert spaces of the subsystems. 

Having split each single-particle state into an A-part 
and a B-part allows to construct the partial densities of 
states (pDOSs) for subsystems A and B and for the two 
spin channels. Summing up the resulting four pDOSs 
over all occupied single-particle states i yields then the 
four electron numbers 

n m 

N A = EEEE/^^ (8a) 

k i j=l k—1 

n n 

n b = EEE E f?<#Sjk<& . (8b) 

k i 3 — 1 k=m+l 

where CT is the occupation number of the single-particle 
Kohn-Sham state i, typically chosen to be a Fermi func- 
tion. With this, also the total spin S of each of the two 
subsystems is determined through the difference of the 
corresponding electron numbers 

S A - \N{-Ni\ (9a) 
Sb = Wl-Ni\ . (9b) 

B. Constraining the Electron Numbers 

Having established the various electron numbers, we 
proceed to introduce the constraint. Aiming to compute 
the non-adiabatic PES . T , ({R/crj}) represent- 

ing a defined spin and charge state in subsystems A and 
B, we first distribute the corresponding iV^, N^, N^, 

and iVg electrons into the four pDOSs derived from the 
set of Kohn-Sham single-particle wavefunctions. This de- 
fines four Fermi energies Cp A , e F A , Cp B and e F B , as well 
as four partial electron densities, which sum to the total 
electron density. If the Fermi energies were all degenerate 
at this stage, the chosen non-adiabatic charge and spin 
state of the system would correspond to the adiabatic 
ground state. However, typically the four Fermi energies 
are all different, reflecting the fact that there is no self- 
consistency between the total electron density and the 
effective Kohn-Sham potential. 

In order to achieve this self-consistency, we choose to 
align the Fermi energies, still requiring that the elec- 
tron numbers that can be filled into the four pDOSs up 
to the resulting common Fermi level remain unchanged. 
We therefore employ a method that shifts each of the 
four pDOSs independently. This is particularly impor- 
tant when hybridization between the two subsystems is 
present and a single-particle state i contains non-zero ex- 
pansion coefficients c kcr for basis functions belonging to 
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subsystem A, as well as for basis functions belonging to 
subsystem B. In such cases, simply shifting the entire 
state i as is e.g. done in the ASCF method would not 
change the relative positions of the A and B Fermi en- 
ergies. Instead, it is necessary to act separately on the 
basis functions in both subsystems. Only this gives the 
electronic structure enough flexibility to fully relax under 
the imposed constrained electron numbers, which in the 
end will yield a different expansion of state i in terms of 
A and B basis functions. 

In practice, we first align the Fermi energies separately 
for each spin. Specifically, we request that eg A = e F b = 
eg, and appropriately shift the A and B Fermi energies 
to the Fermi energy eg. This is achieved by adding to 
the Kohn-Sham Hamiltonian auxiliary potentials that act 
differently on the A and B basis functions. This auxil- 
iary potentials consist of a "strength" factor T a and a 
projection operator P k onto the sub-spaces 
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where the summation is done over the m basis functions 
spanning the Hilbert space of subsystem A and the n — m 
basis functions spanning the Hilbert space of subsystem 
B. In the chosen form P k symmetrically contains covari- 
ant and contravariant basis function a 22 ! 23 , x k an d X*' k ! 
respectively. The resulting matrix representations of P k 
and Pg in the Hilbert space spanned by the Bloch basis 
functions is then hermitian, which facilitates the imple- 
mentation into existing DFT codes as further described 
below. As derived in the Appendix for subsystem A, the 
form of these matrices is as shown schematically in Fig. 
[TJ The matrix element P k ■ • is equal to the overlap ma- 
trix element Sh, if i and j both refer to basis functions 
assigned to subsystem A. If only i or j refer to a basis 
function assigned to subsystem A, the matrix element is 
| • Sfj. Finally, if neither i nor j belong to subsystem 
A, P k j ■ is zero, reflecting that the pDOS of subsystem 
B is not affected by the auxiliary potential r^P k . The 
auxiliary potential rg,P k has an analogous structure. 

Since the purpose of the auxiliary potentials is to align 
the Fermi energies of subsystems A and B with eg, the 
obvious choice for the "strength" factors and Tg is 
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However, because of the resulting non-zero auxiliary po- 
tentials, the initial single-particle states are no longer 
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FIG. 1: Matrix representation of the projection operators P\ 
and Pb in the Hilbert space of the Bloch basis functions %j ■ 
The Sfj are the overlap matrix elements between basis func- 
tions i and j. n is the total number of basis functions and m 
the number of basis functions of subsystem A. 



solutions of the new effective Hamiltonian for each spin, 
comprised of both the Kohn-Sham Hamiltonian and the 
auxiliary potential. As a result, and rg must be de- 
termined self-consistently. Diagonalization of the new ef- 
fective Hamiltonians yields new eigenvectors to construct 
new partial densities of states, the Fermi levels of which 
define new strength factors through Eq. (fTTj) . This is 
repeated in a self-consistency (sc) cycle, until the Fermi 
energies of subsystems A and B are aligned to an arbi- 
trary precision (for each spin channel), 
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The ensuing step of aligning the two different spin 

T I 

Fermi energies e F and e F is done in an analogous way, 
i.e. by adding another auxiliary potential of the form of 
Eq. (jT0|> . In this case, the matrix structure of the corre- 
sponding projection operator P k onto one spin sub-space 
is simpler though. Since this sub-space is spanned by all 
n basis functions, regardless of whether they are in sub- 
system A or B, the sum in Eq. (fTOj) goes up to n, and the 
matrix representation of P k becomes simply the overlap 
matrix, cf. Fig. [T]with m = n. Adding an auxiliary po- 
tential r CT P k to the effective Hamiltonian resulting from 
the preceding alignment of the A and B Fermi energies, 
corresponds therefore to a mere shift of the eigenvalues, 
depending on the chosen "strength" factor T a . Here we 
choose to shift the spin-up and spin-down pDOSs in op- 
posite directions using Ae F = jfi e F~ £ f)' T T = +N^-Aep 
and T l = -A^-Aep with N 1 * = N^+N^, N l = N^ + N^ 
and N = N 1 ' + NK As before, this procedure has to be 
done in a self-consistent way, since adding the new auxil- 
iary potential modifies the effective Hamiltonian. In fact, 
since the alignment of the A and B Fermi levels and the 
alignment of the spin-up and spin-down Fermi levels is 
not independent of each other, the two sc cycles must be 
nested. Discussing below how this can be implemented in 
a numerically efficient way into existing DFT codes, we 
note that once the double self-consistency is achieved, 
we arrive at a final common Fermi level in the system 
and a new set of single-particle states i (with eigenvalues 
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e'f and eigenvectors ip'f) that is the self-consistent solu- 
tion to the effective Hamiltonian, containing the original 
Kohn-Sham Hamiltonian and the two auxiliary poten- 
tials. In each sub-space spanned by the Bloch basis func- 
tions at one k-point we therefore have 

= [h£ s + rx iScf p£ + rg iScf p£ + r* f P k ] ^ 

= f£*$* > (13) 

with "strength" factor values scf , Tg scf and r^ cf as 
determined in the last cycle of the nested sc loops. 



C. Numerical Considerations 

At this stage it is appropriate to recall what has been 
achieved so far. The imposed constraint of fixed elec- 
tron numbers Nl, N^, iVg, and N-q has been suitably 
transformed into external potentials. Adding these to the 
Kohn-Sham Hamiltonian led to the effective Hamiltonian 
of Eq. (fT3|) . As proven by the Hohenberg-Kohn the- 
orem, the electron density resulting from occupying all 
single-particle states up to the common Fermi level cor- 
responds therefore to the fully relaxed electronic struc- 
ture under the given constraint. Calculating the total 
energy, E e T , T , , connected to this electron density 

provides then (together with y nuc ~ nuc ) th e non-adiabatic 
PES U"t „.j »,t „.j i representing the chosen spin and 

electron numbers in the two subsystems. Additionally, 
the derivatives of the non-adiabatic PES with respect to 
the atomic positions can be obtained analytically in the 
standard way providing the forces acting on the atoms. 

At first glance, it appears as if our LC-DFT formalism 
builds on the self-consistent solutions to the Kohn-Sham 
Hamiltonian and then requires two additional nested sc 
loops. However, for the latter self-consistency under the 
imposed constraint there is no need to start with self- 
consistent Kohn-Sham solutions. Additionally, the mere 
eigenvalue shift induced by the spin Fermi level alignment 
step is nicely compatible with the structure of existing 
DFT codes. From a computational point of view, our 
algorithm can thus be implemented as one additional sc 
cycle at each iteration of the existing electronic sc cycle in 
a DFT code. Of course, all known approaches to improve 
the convergence of sc cycles like the use of sophisticated 
mixing schemes can equally be applied to the additional 
cycle. Having implemented a Pulay mixing scheme 24 , our 
experience was that for the tested systems only the very 
first DFT iterations required a significant number of in- 
ner iterations, typically about 5 to 10, while close to the 
outer self-consistency also the inner self-consistency was 
often directly reached after the first iteration. For the ex- 
ample involving the Al(lll) surface, we even found the 
overall DFT convergence, i.e. the number of iterations in 
the outer sc cycle, frequently much improved compared 
to standard adiabatic calculations, because the oscilla- 
tion of states around the Fermi level is reduced when 



controlling the electron numbers and thereby also the 
occupation numbers of the Kohn-Sham states close to 
the Fermi level. For the systems with localized electronic 
states, we found that using the spins or electron num- 
bers in the subsystems as convergence criterion for the 
self-consistency was sometimes preferred to the equality 
of the Fermi energies. Overall, for the systems studied, 
the cost of the calculations using the constraint was thus 
often about the same and sometimes even lower than the 
cost of standard adiabatic calculations. 



D. Comparison to ASCF and Fixed-Spin Moment 

An important characteristic of the LC-DFT approach 
presented here is that it is parameter-free and only the 
set of the electron numbers in the four channels defining 
the non-adiabatic state of interest has to be specified. 
Under this constraint, the electronic structure is fully re- 
laxed, i.e. the partial densities of states are not frozen 
and shifted statically. Instead, the single-particle states 
are flexible to vary the contribution of each basis func- 
tion to each single-particle state freely. This can lead to a 
significant improvement compared to the two prominent 
and widely employed implementations of the constrained 
DFT concept, the ASCFii and the fixed -spin moment 
approach 1 -. In the latter approach the system is only 
separated into a spin-up and a spin-down channel, which 
are filled independently. In LC-DFT we go a step further 
and allow also to distribute the spin-up and spin-down 
electrons in a well-defined way into the two subsystems, 
which permits a more general control over the spatial 
distribution of the electron and magnetization densities. 
The different results obtained with both methods are par- 
ticularly obvious for the oxygen dissociation at Al(lll) 
case described in Section IIIC below. 

The improvement compared to the ASCF method con- 
cerns extended systems. Both approaches have in com- 
mon that they consider two subsystems and first ana- 
lyze the pDOSs to determine the contributions of each 
subsystem to the individual single-particle states. How- 
ever, in the ASCF method the states with high A-parts 
are then completely assigned to subsystem A, regardless 
of their B-contributions, while the states with small A- 
contributions are completely assigned to subsystem B^& 
A subsequent occupation of the states by or 1 elec- 
trons is therefore only fully justified for the typically not 
interesting case of non-interacting subsystems, i.e. no 
hybridization. Otherwise it results in fractional effec- 
tive occupation numbers of the individual subsystems, 
which necessarily introduces some uncertainty in the to- 
tal energies obtained in the ASCF method. The LC-DFT 
approach, on the other hand, allows for a physical rehy- 
bridization of the states under the imposed constraint It 
thus conserves the electron numbers of both subsystems, 
while at the same time fully taking hybridization into ac- 
count. This is the main difference between the present 
LC-DFT formalism and the ASCF method, while in the 
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limit of infinitely separated subsystems both approaches 
are equivalent. 

It is finally also important to note that all methods 
discussed here, the FSM approach, the ASCF method 
and the LC-DFT formalism, intend to overcome limita- 
tions in the description of chemical processes by control- 
ling the electron numbers. This does obviously not allow 
to overcome approximations in the employed exchange- 
correlation (xc) functional. Local-density or gradient- 
corrected xc functionals are e.g. known to cause inac- 
curacies in the energy splittings between different spin 
multiplet o 27 i 28 ' 29 . To this end, we note that in the limit 
of infinite separation between the subsystems (reactants), 
our approach reduces to the description of two isolated 
subsystems. The non-adiabatic states of defined spin or 
charge in our approach are then simply the correspond- 
ing excited states of the isolated subsystems, i.e. of the 
non-interacting reactants. It is worthwhile to point out 
that this is different to the diabaticity concept, which 
reduces in this limit of infinite separation to the adia- 
batic (excited) states of the interacting system. If there 
is a (unphysical) charge or spin transfer even at these 
distances (as in the example of 02/Al(lll) discussed be- 
low), it will be present in both the adiabatic and the dia- 
batic description, but not in our non-adiabatic approach. 
The reduction to the states of the non-interacting reac- 
tants provides also a possibility to check the accuracy and 
suitability of the employed approximate xc functional 
to describe the system states of interest. For the non- 
interacting reactants it is typically feasible to check the 
results obtained with the employed xc functional against 
higher- level calculations. If one finds a sufficiently accu- 
rate description there, e.g. for the lowest state of differ- 
ent symmetries, then we believe it is a valuable approach 
to maintain this functional also for closer distances of 
the reactants and get approximate insight into the cor- 
responding non-adiabatic potential-energy surface. 



III. APPLICATIONS 

As an important application of the LC-DFT method 
we now illustrate its use in the calculation of non- 
adiabatic PESs that are of interest in the investigation 
of dynamic processes. Specifically, we focus here on the 
scattering of atoms or molecules in crossed molecular 
beams or at solid surfaces. As a side effect this also 
shows how the method can be employed to suitably re- 
strict unwanted electron transfer between weakly or non- 
interacting subsystems. In a DFT calculation, such an 
electron transfer will e.g. occur whenever an occupied 
level of subsystem A is higher in energy than an unoccu- 
pied level of subsystem B. Alignment of the Fermi levels 
of the two systems will then lead to a fractional occupa- 
tion of both states in the self-consistent solution, even 
if the distance between the two subsystems is macro- 
scopic. A prominent example for this is the interaction 
of an oxygen molecule with the Al(lll) surface, where 



Na + CI 




JV Na JV Na 






Ionic 


5 5 


9 


9 


Neutral, singlet 


6 5 


8 


9 


Neutral, triplet 


6 5 


9 


8 


Mg 4 + Cl 2 




/V 1 N L 

JV M g4 JV M g4 






Neutral 


24 24 


17 


17 


Ionic, singlet 


23 24 


18 


17 


Ionic, triplet 


24 23 


18 


17 


2 + Al(lll) 




Nl Nk 


JV A1(111) 


N L 

ly Al(lll) 


Neutral, triplet 


18 14 


409.5 


409.5 



TABLE I: Constrained electron numbers for the three test 
systems discussed. 

the unoccupied 27r*^ orbitals of the molecule are lower 
than the Fermi level of the metah 17 i 25 The resulting elec- 
tron transfer lowers the total energy of the system and at 
macroscopic distances the latter does then not converge 
to the physically correct limit, given by the sum of the 
total energies of the isolated molecule and the isolated 
surface. 

As an example for an extended periodic system, we cor- 
respondingly briefly discuss the interaction of O2 with the 
Al(lll) surface. In addition, we also present calculated 
non-adiabatic PESs for two finite model systems, namely 
for the scattering of Na and CI, and the scattering of CI2 
at a small Mg4 cluster. All calculations have been carried 
out using the all-electron DFT code DMol 3 , which em- 
ploys numerical atomic-like orbitals as basis functions* 3 ^ 
Unless otherwise noted, we employ the DMol 3 standard 
basis set "all" consisting of atomic orbitals and polariza- 
tion functions 3 ^, a real-space cutoff of 12 bohr for the 
basis functions, and the PBE 3 -^ xc functional. 

A. Na + CI 

In the scattering of a sodium and a chlorine atom, two 
non-adiabatic states of interest are the ionic PES "Na + 
+ CI"" and the neutral PES "Na + CI". Since both neu- 
tral atoms are spin doublets in their ground states, there 
are two possible relative orientations of the spins in the 
latter case, yielding an overall singlet (antiparallcl spins) 
or triplet (parallel spins) neutral state. Identifying each 
atom as one subsystem in our LC-DFT approach, Table 
U shows the constrained electron and spin numbers we 
used to represent each of these non-adiabatic states. The 
resulting PESs as a function of the interatomic separa- 
tion are shown in Fig. [21 in which we additionally include 
the computed adiabatic ground state PES and the PES 
as resulting from a FSM calculation for an overall spin- 
triplet state (iV T = 15, N l = 13). In both the latter 
PES and the neutral triplet PES there are therefore two 
unpaired electrons, but only the neutral triplet PES com- 



7 



-2 




2 4 6 8 10 

Z(A) 

FIG. 2: (Color online) Non-adiabatic potential-energy sur- 
faces (PESs) for the scattering of a Na and a CI atom. Shown 
are the energies as a function of the interatomic distance Z. 
See Table |T] and text for the constrained electron numbers 
defining the various non-adiabatic PESs. Additionally shown 
are the adiabatic ground state PES and the PES obtained 
with a fixed-spin moment (FSM) triplet calculation. The en- 
ergy zero corresponds to the energy of the two isolated, neu- 
tral atoms. 



puted by LC-DFT has the additional control of locating 
one unpaired electron explicitly at each atom. 

For bond lengths lower than 8 A, the energetically most 
favorable non-adiabatic state is found to be the ionic 
PES, whereas for larger distances these are the degen- 
erate singlet and triplet neutral PESs. The degeneracy 
of the latter two PESs is only lifted for distances smaller 
than 5 A, at which the Pauli repulsion between the two 
unpaired spin-up electrons leads to a strong increase of 
the neutral triplet PES. Interestingly, the FSM curve 
is for all distances virtually degenerate to this neutral 
triplet PES, indicating that even without constraint one 
unpaired electron wants to stay at each atom. The min- 
imum of the adiabatic PES, on the other hand, is some- 
what lower than the minimum of the ionic PES, showing 
that the electron transfer in the adiabatic case differs 
slightly from the one electron imposed in the LC-DFT 
computation. 

Another prominent feature of Fig. [2] is that for large 
interatomic separations the adiabatic PES is about 1 eV 
lower than the limit of neutral separated atoms. This 
is an illustration of the above described electron transfer 
problem in adiabatic calculations. Even for infinite inter- 
atomic distances, a small amount of charge is transferred 
from the sodium to the chlorine atom, since the lowest 
unoccupied state (LUMO) of the latter, is lower in 
energy than the highest occupied 3s^ state (HOMO) of 
the sodium atom. In the self-consistent calculation, elec- 
tron density is consequently transferred, until the Fermi 
levels of the two atoms are aligned. At infinite separation 
between the two atoms, this charge transfer can be deter- 
mined quantitatively by calculating the Na Kohn-Sham 
HOMO and CI Kohn-Sham LUMO level energies as a 
function of different occupations. The results obtained 




0.1 0.2 0.3 0.4 0.5 
charge transfer from Na to CI (e) 



FIG. 3: (Color online) Level energy of the highest occupied 
Kohn-Sham molecular orbital (HOMO) of a free Na atom 
and of the lowest unoccupied Kohn-Sham molecular orbital 
(LUMO) of a CI atom as a function of the electronic occupa- 
tion. The occupations in the two separate calculations of the 
isolated atoms are varied in form of a charge transfer, i.e. the 
Na atom is computed with a fraction of an electron removed, 
and the CI atom is computed with this fraction of an electron 
added. 




2 4 6 8 10 

z(A) 

FIG. 4: Force acting on the CI atom in the NaCl dimer as 
a function of the interatomic distance Z for the ionic PES. 
Filled squares represent the forces calculated within the LC- 
DFT approach, open triangles the forces resulting from a nu- 
merical differentiation of the PES shown in Fig. [2] 



from calculations of the isolated charged atoms are dis- 
played in Fig. [3] and show that HOMO and LUMO are 
only aligned after an electron transfer of 0.37 e. This 
unphysical electron transfer at infinite separations is not 
possible in the LC-DFT approach by construction, ex- 
plaining why in contrast to the adiabatic PESs the neu- 
tral triplet and singlet PESs approach the correct limit. 

Finally, we also employed this system to test the proper 
evaluation of the forces in the constrained LC-DFT cal- 
culations. Taking the example of the ionic PES, Fig. 0] 
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FIG. 5: (Color online) Non-adiabatic potential-energy sur- 
faces (PESs) of a CI2 molecule scattering at a Mg4 cluster. 
Shown are the energies as a function of the distance Z of the 
centers-of-mass of the CI2 molecule and the Mg4 cluster, with 
a scattering geometry as explained in the inset. See Table 
U and text for the constrained electron numbers defining the 
various non-adiabatic PESs. Additionally shown is the adi- 
abatic ground state PES, and the energy zero corresponds 
to the energy of the isolated neutral CI2 molecule and Mg4 
cluster. 



shows the force on the CI atom computed either ana- 
lytically within the LC-DFT approach or numerically by 
differentiating the PES shown in Fig. O The agreement 
is excellent proving that forces can be computed accu- 
rately within the LC-DFT approach. 



B. Mg 4 + Cl 2 

As a simple example for subsystems consisting of 
groups of atoms we discuss the interaction of a CI2 
molecule with a small metal cluster formed of a tetra- 
hedron of 4 magnesium atoms. To illustrate the method 
we restrict ourselves here to computing the PESs as a 
function of the distance of a CI2 molecule approaching 
the Mg3 plane of the cluster as explained in Fig. [5] For 
this we first relaxed the structure of the cluster and the 
CI2 molecule separately, and then held the resulting Mg- 
Mg bond lengths of 3.07 A and the Cl-Cl bond length 
of 2.03 A fixed in the subsequent calculations. Using 
the electron numbers compiled in Table [H we calculated 
the non-adiabatic PESs corresponding to the neutral, the 
ionic-singlet, and the ionic-triplet state. Together with 
the adiabatic PES, the resulting curves are shown in Fig. 
[5j At all distances, the non-adiabatic state corresponding 
to a neutral configuration exhibits the lowest energy, with 
the two ionic curves exhibiting significantly higher ener- 
gies. Similar to the Na+Cl case the latter two become 
degenerate at larger distances, when the Pauli repulsion 
affecting the ionic triplet curve becomes negligible. The 
closeness of the non-adiabatic neutral curve to the adi- 
abatic result indicates only a comparably small electron 



transfer from the cluster to the molecule during the scat- 
tering process. In this case, there is therefore also only a 
small electron transfer problem and the adiabatic curve 
approaches the proper limit for large molecule-cluster 
separations. By differently occupying the Kohn-Sham 
HOMO and LUMO levels as done in Fig. [31 we indeed 
obtain only a very small electron transfer of 0.02 e that 
is required to align the Fermi energies in this case. 



C. O2 Dissociation at Al(lll) 

As a final example for an extended system, treated by 
periodic boundary conditions, we turn to the dissocia- 
tion of an O 2 molecule at the Al(lll) surface. For this 
system the postulated dominant role of non-adiabatic 
effect o 17 ' 25 i 32 could not be verified until recently, since 
only empirical estimates of the underlying non-adiabatic 
PESs were available^ 3 - 3 ^. Using LC-DFT we now fo- 
cus on the non-adiabatic neutral triplet PES, which is a 
suitable representation at large distances from the sur- 
face, where the gas-phase O2 molecule will be in its spin- 
triplet ground state and the Al(lll) surface in a spin- 
singlet state. For the calculations we employed a (3x3) 
Al(lll) slab consisting of 7 aluminium layers, separated 
by a 30 A vacuum. Oxygen is adsorbed at both sides 
of the slab to establish inversion symmetry, a real space 
cutoff of 9 bohr has been applied to the basis functions, 
and 10 k-points have been used to sample the irreducible 
wedge of the Brillouin zone. The electron numbers of the 
oxygen and the aluminium subsystem used to define the 
neutral triplet PES are listed in Table |TJ 

Discussing our results for the high-dimensional PES 
in detail elsewher e) 17 ! 36 , we illustrate the insights gained 
by the LC-DFT approach by concentrating on the two- 
dimensional dependence on the molecular bond length 
r and the center-of-mass distance of the molecule from 
the surface Z for a fixed molecular orientation and lateral 
position over the surface. Fig. [J5] shows corresponding "el- 
bow plots" specifically for an O2 molecule approaching 
the surface head-on and above an fee site. In agreement 
with previous studie a 25 ' 32 the adiabatic PES displayed 
in Fig. [6^ does not exhibit an energy barrier to dissoci- 
ation, a finding that cannot be reconciled with the ex- 
perimentally well-established low sticking coefficient for 
thermal molecule s 17 ' 37 . Suspecting a dominant role of 
non-adiabatic effects as the reason for this discrepancy, 
we turn to non-adiabatic representations, in which par- 
ticularly the spin-triplet character of the gas-phase O2 
molecule is conserved. Fig. and c show correspond- 
ing PESs obtained with the FSM approach for an overall 
triplet state of the system and with our LC-DFT ap- 
proach constraining the spin-triplet to the O2 molecule, 
respectively. In contrast to the Na+Cl neutral triplet 
PES discussed above, the FSM and LC-DFT approach 
now yield qualitatively different results. While no barrier 
is obtained in the prior method, the neutral triplet PES 
calculated with LC-DFT exhibits a clear energy barrier. 
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r (A) r (A) r (A) 

FIG. 6: (Color online) Two-dimensional cut ("elbow plot") through the high-dimensional PES for the O2 dissociation at 
Al(lll). The energy is shown as a function of the center-of-mass distance of the molecule from the surface Z and the oxygen- 
oxygen bond length r. The molecule approaches the surface head-on above an fee site as explained in the insets, (a) Adiabatic 
calculation, (b) triplet fixed-spin moment calculation, (c) neutral triplet LC-DFT calculation. Only the latter PES exhibits 
an energy barrier. The energy difference between the contour lines is 0.2 eV and the small white circle denotes the molecular 
position discussed in Fig. 7. 




FIG. 7: (Color online) Magnetization density (difference be- 
tween spin-up and spin-down electron density) for a molecule 
at the energy barrier of the LC-DFT triplet PES (r = 1.4 A, 
Z = 1.8 A, as marked in Fig. [6}. Shown is a two-dimensional 
cut perpendicular to the surface and through the O2 molecule, 
along the [Oil] and [111] direction. The position of the two 
O atoms are marked as small black circles, the position of 
the Al atoms as small white circles, (a) Isolated O2 molecule 
in its spin-triplet ground state, (b) adiabatic calculation for 
Oa/Al(lll), (c) triplet FSM calculation for 2 /Al(lll), and 
(d) neutral triplet LC-DFT calculation for 2 /Al(lll). 



the O2 molecule, nor the metal atoms exhibit any spin- 
polarization, which is the most favored state for small 
molecule-surface separations. In the FSM calculation, 
the total spin of the system is forced to be a triplet, but 
as apparent from (c) the majority of this excess spin is 
not located at the O2 molecule, but distributed over the 
entire metal slab. In contrast to this, in the LC-DFT re- 
sult shown in (d) the triplet spin is localized at the oxygen 
molecule, reflecting the improved control over the spatial 
distribution of the magnetization density in the latter 
approach. The accumulation of spin-up density on the 
O2 molecule repels the spin-up density of the metal slab 
towards the interior of the slab, while there is a strong 
accumulation of spin-down density close to the molecule. 
As a consequence, the metal slab is still in an overall 
singlet state in the LC-DFT calculation, which is obvi- 
ously a better representation of the non-adiabatic state 
defined by an impinging triplet O2 molecule compared to 
the FSM result. In addition, the LC-DFT approach over- 
comes the small charge transfer problem present in this 
system, as well. Since the unoccupied 2tt*^ orbitals of the 
O2 molecule are lower than the Fermi level of the metal, 
the adiabatic calculation yields a charge transfer of 0.01 e 
to the O2 molecule even at macroscopic distances from 
the surface. 



The reason for this difference is the different distribu- 
tion of the magnetization density in the system. This 
is illustrated for a molecular configuration at the en- 
ergy barrier in Fig. [7J In (a) the magnetization den- 
sity computed for the free oxygen molecule (i.e. without 
Al(lll) slab) in its spin-triplet ground state is shown, 
whereas (b) displays the result of an adiabatic calcula- 
tion including the Al(lll) slab. In the latter case, neither 



IV. SUMMARY 

In summary, we presented a locally-constrained 
density-functional theory (LC-DFT) approach, that al- 
lows to confine electrons to sub-spaces of the Hilbert 
space, e.g. to selected atoms or groups of atoms. A 
major application of this technique is the computation of 
non-adiabatic potential-energy surfaces, which we illus- 
trated with examples treating the scattering of atoms and 
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molecules at other atoms, clusters or surfaces. Following 
the general formulation by Dederichs et alJ^, the elec- 
tron confinement is realized by suitably introducing ad- 
ditional constraints to the electronic Hamiltonian. DFT 
is then used to obtain the fully relaxed electronic struc- 
ture under the additional external potential imposed by 
the applied constraints. With the ASCF and FSM meth- 
ods as widely employed alternative implementations of 
this general concept, our LC-DFT method offers a more 
systematic approach to extended systems compared to 
ASCF, and better control over the spatial distribution of 
the constraint electrons compared to FSM. This better 
spatial control allows also to overcome the charge trans- 
fer problem between widely separated subsystems that 
can occur in adiabatic DFT calculations. 
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Appendix: The Projection Operator 

In general, localized atom-centered basis functions like 
atomic orbitals are not orthogonal, and a projection op- 
erator should be formulated in terms of covariant and 
contravariant basis functions^.^. The covariant Bloch 
basis functions \xf > and the contravariant Bloch basis 
functions |x jk > are related to each other by the equa- 
tions 



alent, 



< x 3k \ 
\x 3k > 
<x k l 

and \Xi > 



k\-l 



1=1 

n 

E Ix? > 

1=1 

n 

E sk <* ik i > 



3=1 



E i^' k > ^ . 

3=1 



(14a) 
(14b) 
(14c) 
(14d) 



with S k being the overlap matrix. For covariant and 
contravariant basis functions the following orthonormal- 
ity relation holds, 



< X*lx k > = < Xi\x jk > = S, 



(15) 



where 5^ is the Kronecker symbol. 

In principle there are two possible forms of the projec- 
tion operator P k into the subsystem A, which are equiv- 



pk 



Pa = Elx*><X k | 



(16a) 



(16b) 



and where the sum i = 1 . . . m runs over the m basis func- 
tions of subsystem A. However, expanded onto the Bloch 
basis functions, both these formulations for P k yield non- 
hermitian matrices. In order to facilitate the implemen- 
tation into existing DFT codes, we therefore prefer to 
work with the following symmetrized form, which does 
lead to a hermitian matrix 



pk 

*A 



§• Elx k xx ik l + Elx 

\i=l i=l 



k >< x k l 



Inserting the expressions for the contravariant basis 
functions in Eq. (| 14[) enables us to express the projection 
operator entirely in terms of the known covariant basis 
functions (the atomic orbitals) 



1 



m n 



(18) 



»=i j=i 

m n 



+ EEix3 k >( 5k i) _1 <x k i 
»=i j=i 

Starting from this expression, we can then derive the ma- 
trix representation of this operator in the Hilbert space 
spanned by the Bloch basis functions 

< x k |P k lx k > 



pk 

^A,ij 



/ rn n 

= 2<X k | EElx k >(^ k ) _1 <X k | 

\fe=l 1=1 

tn n \ 

+ EEix k >( sk r 1 <x k i ix k > 

r=l s=l / 
1 / m n 

= 2 EE< ^ > < x k ix k > 
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m n \ 

+ EE<A"> ofc) -1 < x k ix k > 
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_^ / m n 

= ^(EE^^^^ 1 ^ k 

Vfe=i 1=1 
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+ &is i^sr ) ^Srj J 

r=l s=l / 
^ / m m \ 

= 2 (E s,k<5fe j + E^ 5 ^' J 



{Sfj for i < m A j < m 
T^Sfj for i < m V j < m 
for i > m A j > m 



(19) 
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